BarnesOI Subroutine

SUBROUTINE BarnesOI &
!
(network, grid, gammapar)

USE Units, ONLY: &
!Imported parameters
pi

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network

REAL (KIND = float), OPTIONAL, INTENT(IN) :: gammapar

!Arguments with intent (inout):
TYPE (grid_real), INTENT (INOUT) :: grid

!local declarations:
REAL (KIND = float) :: dn         !!data spacing used in the analysis
REAL (KIND = float) :: gamma      !! smoothing parameter [0.2 - 1.0]
REAL (KIND = float) :: xkappa_1   !! Gauss constant for first pass
REAL (KIND = float) :: xkappa_2   !! Gauss constant for second pass
REAL (KIND = float) :: rmax_1     !! maximum scanning radii, for first
REAL (KIND = float) :: rmax_2     !! and second passes
REAL (KIND = float) :: anum_1     !! numerator, beyond scanning radius,
REAL (KIND = float) :: anum_2     !! for first and second passes
REAL (KIND = float) :: xa,ya      !! x, y coords of current station
REAL (KIND = float) :: xb,yb      !! x, y coords of current station
REAL (KIND = float) :: w1,w2      !! weights for Gauss-weighted average
REAL (KIND = float) :: wtot1,wtot2 !!sum of weights
REAL (KIND = float) :: ftot1,ftot2 !!accumulators for values, corrections
REAL (KIND = float) :: dsq !!square of distance between two points
REAL (KIND = float) :: dnMax, dnMin
REAL (KIND = float) :: nc, nr !! number of columns and rows in the grid
REAL (KIND = float) :: cellsize ![m]
REAL (KIND = float), ALLOCATABLE :: dvar (:) !!estimated error
TYPE (ObservationalNetwork) :: activeNetwork
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: i, j, k !! loop indexes

!----------------end of declarations-------------------------------------------

!check for available measurements
CALL ActualObservations (network, count, activeNetwork)

!Allocate dvar
ALLOCATE (dvar(activeNetwork % countObs))

!set points
point1 % system = grid % grid_mapping
point2 % system = grid % grid_mapping

!set gamma to default or user specified value
IF (PRESENT(gammapar)) THEN
  gamma = gammapar
ELSE
  gamma = 0.2 !default value
END IF

!get dn
!compute average cellsize.
cellsize = ( CellArea (grid, grid % idim / 2, grid % jdim / 2) ) ** 0.5 ![m]
nc = grid % jdim
nr = grid % idim
dnMax = ( (nr * cellsize) * (nc * cellsize) ) ** 0.5 * &
        ((1.0 + activeNetwork % countObs ** 0.5) / (activeNetwork % countObs - 1.0))

dnMin = ( (nr * cellsize) * (nc * cellsize) / activeNetwork % countObs ) ** 0.5

dn = 0.5 * (dnMin + dnMax)
!debug
dn = 6. * cellsize

! Compute the first and second pass values of the scaling parameter
! and the maximum scanning radius used by the Barnes scheme.
! Values above this maximum will use a 1/r**2 weighting.

! First-round values
	xkappa_1 = 5.052 * (2. * dn / Pi) ** 2.0

! Define the maximum scanning radius to have weight defined by
! wt = 1.0 x 10**(-30) = exp(-rmax_1/xkappa_1)
! Also scale the 1/r**2 wts so that when dsq = rmax, the wts match.
	rmax_1 = xkappa_1 * 30.0 * log(10.0)
	anum_1 = 1.E-30 * rmax_1
	
! Second-round values
	xkappa_2 = gamma * xkappa_1
	rmax_2 = rmax_1 * gamma
	anum_2 = 1.E-30 * rmax_2
	
! Scan each input data point and construct estimated error, dvar, at that point
outer_loop: DO i = 1, activeNetwork % countObs

        point1 % northing = activeNetwork % obs (i) % xyz % northing  
        point1 % easting  = activeNetwork % obs (i) % xyz % easting
		wtot1 = 0.0
		ftot1 = 0.0
		
   inner_loop: DO j = 1, activeNetwork % countObs
   
            point2 % northing = activeNetwork % obs (j) % xyz % northing  
            point2 % easting  = activeNetwork % obs (j) % xyz % easting
            dsq = Distance (point1,point2) ** 2.
            
            IF (dsq <= rmax_1) THEN
               w1 = exp((- dsq)/xkappa_1)
            ELSE !Assume a 1/r**2 weight
               w1 = anum_1/dsq
            END IF

			wtot1 = wtot1 + w1
			ftot1 = ftot1 + w1 * activeNetwork % obs (j) % value
    END DO inner_loop

!
!		if(wtot1==0.0) printf("stn wt totals zero\n");

        dvar(i) = activeNetwork % obs (i) % value - ftot1/wtot1

END DO outer_loop

! Grid-prediction loop.  Generate the estimate using first set of
! weights, and correct using error estimates, dvar, and second
! set of weights.

row_loop: DO i = 1, grid % idim
  col_loop: DO j = 1, grid % jdim
  
  IF (grid % mat (i,j) /= grid % nodata) THEN
      CALL GetXY (i,j,grid, point1 % easting, point1 % northing)

      ! Scan each input data point.
      ftot1 = 0.0
      wtot1 = 0.0
      ftot2 = 0.0
      wtot2 = 0.0
      
      stations_loop: DO k = 1, activeNetwork % countObs
          point2 % northing = activeNetwork % obs (k) % xyz % northing  
          point2 % easting  = activeNetwork % obs (k) % xyz % easting
          
          dsq = Distance (point1,point2) ** 2.

          IF (dsq<=rmax_2) THEN
              w1 = exp((- dsq)/xkappa_1)
              w2 = exp((- dsq)/xkappa_2)
          ELSE IF (dsq<=rmax_1) THEN
              w1 = exp((- dsq)/xkappa_1)
              w2 = anum_2/dsq;
          ELSE  
              !Assume a 1/r**2 weight.
              w1 = anum_1/dsq
              !With anum_2/dsq.
              w2 = gamma * w1
          END IF

          wtot1 = wtot1 + w1
          wtot2 = wtot2 + w2
          ftot1 = ftot1 + w1 * activeNetwork % obs (k) % value
          ftot2 = ftot2 + w2 * dvar (k)
 
      END DO stations_loop
!
!			if (wtot1==0.0 || wtot2==0.0) printf("wts total zero\n");
!
			grid % mat (i,j) =  ftot1/wtot1 + ftot2/wtot2

  END IF
  END DO col_loop
END DO row_loop



DEALLOCATE(activeNetwork % obs)
DEALLOCATE (dvar)


RETURN
END SUBROUTINE BarnesOI